program main
! this is a fokker-planck code  solving spitzer-harm equation to calculate electrical conductivity.
  use precision,only:p_
  use constants,only: m,pi   ! m is the number of velocity grid number.
  implicit none
  integer:: i,j,nstep
  real(p_):: dv,dt,zi
  real(p_):: v(m),q0(m),q(m),fm(m)
  real(p_):: aa(m),bb(m),cc(m),col(m),dvv(m),djj(m)
  real(p_):: a(m),b(m),c(m),r(m)
  real(p_):: cond     !electrical conductivity
! this namelist provides the parameter Zi and the step number in time.
  namelist /input/ nstep,zi  
  open(40,file='input_file')
  read(40,input)
  !create grid and give maxwellian distribution and initial distribution for q.
  dv=12.0_p_/(m-1)
  do j=1,m
     v(j)=0.0_p_+(j-1)*dv
     fm(j)=1./sqrt((2.*pi)**3)*exp(-0.5*v(j)**2)
     q0(j)=1.0 
  enddo
  !time step size, here is set to 10 collision time. 
  !Larger size can be used thanks for the implicit scheme.
  dt=10.0_p_   
  q=q0 ! initial condition for q
  open(30,file='abc.dat')  ! this file stores the steady distribution function.

  ! advance in time
  do i=1,nstep
     !calculate the difusion coefficient for maxwellian background.
     call get_d_f(dv,v,fm,dvv,djj) 
     ! calculate the coefficient of the differential equation (adjoint equation)
     call get_coefficient(zi,dv,v,fm,dvv,djj,aa,bb,cc) 
     ! calculate the collision term between maxwellian and non-maxwelian background.
     call get_col(dv,v,fm,q,col)
     ! constuct the tridiagonal matrix
     call get_matrix(dt,dv,v,q,aa,bb,cc,col,a,b,c,r)
     ! sovle the tridiagonal matrix equation
     call tridag(a,b,c,r,q,m)
  enddo
  !output the steady distribution function
  do j=1,m
     write(30,*) v(j),q(j),q(j)*fm(j)
  enddo
  close(30)
  ! calculate the electrical conductivity
  call get_conductivity(dv,v,fm,q,cond)
  write(*,*) 'conductivity=', cond
end program main


subroutine get_d_f(dv,v,fb,dvv,djj)
  ! This routine calculates the difusion coefficient for maxwellian background.
  use precision,only:p_
  use constants,only: m,pi
  implicit none
  real(p_),intent(in):: dv,v(m),fb(m)
  real(p_),intent(out)::dvv(m),djj(m)
  integer::i,j
  real(p_):: sum,sum1,sum2

  do i=2,m-1
     sum1=0.
     do j=2,i
        sum1=sum1+v(j)**4*fb(j)*dv
     enddo
     sum2=0.
     do j=i+1,m
        sum2=sum2+v(j)*fb(j)*dv
     enddo
     dvv(i)=4.*pi/3.*(sum1/v(i)**3+sum2)
  enddo

  do i=2,m-1
     sum1=0.
     do j=2,i
        sum1=sum1+v(j)**2/(2.*v(i)**3)*(3.*v(i)**2-v(j)**2)*fb(j)*dv
     enddo
     sum2=0.
     do j=i+1,m
        sum2=sum2+v(j)*fb(j)*dv
     enddo
     djj(i)=4.*pi/3*(sum1+sum2)
  enddo

end subroutine get_d_f

subroutine get_coefficient(zi,dv,v,fm,dvv,djj,aa,bb,cc)
  !This routine calculate the coefficient of the adjoint equation.
  !aa(m) for coefficient before second derivative
  !bb(m) for coefficient before first derivative
  !cc(m) for coefficient before zero derivative
  use precision,only:p_
  use constants,only: m,pi
  implicit none
  real(p_),intent(in):: zi,dv,v(m),fm(m),dvv(m),djj(m)
  real(p_),intent(out)::aa(m),bb(m),cc(m)
  integer::i,j
  real(p_):: sum1,sum2,tmp

  do i=2,m-1
     aa(i)=dvv(i)
  enddo
  do i=2,m-1
     sum1=0.
     do j=2,i
        sum1=sum1+v(j)**4*fm(j)*dv
     enddo
     sum2=0.
     do j=i+1,m
        sum2=sum2+v(j)*fm(j)*dv
     enddo
     tmp=4*pi/3.*(-1./v(i)**4*sum1+2./v(i)*sum2) 
     bb(i)=tmp-v(i)*aa(i)
  enddo

  do i=2,m-1
     cc(i)=-2.*djj(i)/v(i)**2-zi/v(i)**3
  enddo

end subroutine get_coefficient

subroutine get_col(dv,v,fm,q,col)
  !This routine calculates the collision term between maxwellian and non-maxwellian.
  !This is an implement of Eq.(34) in Karney1986 paper.
  use precision,only:p_
  use constants,only: m,pi
  implicit none
  real(p_),intent(in):: v(m),fm(m),q(m),dv
  real(p_),intent(out)::col(m)
  integer::i,j
  real(p_):: sum1,sum2,tmp

  do i=2,m-1
     sum1=0.
     do j=2,i
        tmp=v(j)**3/(5.*v(i)**2)-v(j)/(3.*v(i)**2)
        sum1=sum1+v(j)**2*tmp*fm(j)*q(j)*dv
     enddo
     sum2=0.
     do j=i+1,m
        tmp=v(i)**3/(5.*v(j)**2)-v(i)/(3.*v(j)**2)
        sum2=sum2+v(j)**2*tmp*fm(j)*q(j)*dv
     enddo
     col(i)=4.*pi*(fm(i)*q(i)+sum1+sum2)
  enddo
end subroutine get_col

subroutine get_matrix(dt,dv,v,q,aa,bb,cc,col,a,b,c,r)
  !This routine is to construct the tridiagonal matrix.
  use precision,only:p_
  use constants,only: m
  implicit none
  real(p_),intent(in):: dt,dv,v(m),q(m),aa(m),bb(m),cc(m),col(m)
  real(p_),intent(out):: a(m),b(m),c(m),r(m)
  integer:: i

  ! the following loop is to calculate 3 diagonal lines in the tridiagonal matrix
  do i=2,m-1
     a(i)=-aa(i)*dt/dv**2+bb(i)*dt/(2.*dv)
     b(i)=1.+2*aa(i)*dt/dv**2-cc(i)*dt
     c(i)=-aa(i)*dt/dv**2-bb(i)*dt/(2.*dv)
     r(i)=q(i)+dt*col(i)+dt*v(i)    !  r(i)=q(i)+dt*col(i)-dt*v(i)

  enddo
  ! boundary conditions:
  ! the following 3 lines are to implement q(1)=0
  b(1)=1._p_
  c(1)=0._p_
  r(1)=0._p_
  ! the following 3 lines are to implement q''(m-1)=0
  a(m)=-b(m-1)-2*a(m-1)
  b(m)=a(m-1)-c(m-1)
  r(m)=-r(m-1)
end subroutine get_matrix

SUBROUTINE tridag(a,b,c,r,u,n)
  use precision,only:p_
  implicit none
  INTEGER n
  integer,PARAMETER:: NMAX=1000
  REAL(p_):: a(n),b(n),c(n),r(n),u(n)
  INTEGER j
  REAL(p_) bet,gam(NMAX)
  if(b(1).eq.0.) stop 'tridag: rewrite equations'
  bet=b(1)
  u(1)=r(1)/bet
  do  j=2,n
     gam(j)=c(j-1)/bet
     bet=b(j)-a(j)*gam(j)
     if(bet.eq.0.) stop 'tridag failed'
     u(j)=(r(j)-a(j)*u(j-1))/bet
  enddo
  do  j=n-1,1,-1
     u(j)=u(j)-gam(j+1)*u(j+1)
  enddo
  !  (C) Copr. 1986-92 Numerical Recipes Software ,4-#.
END SUBROUTINE tridag


subroutine get_conductivity(dv,v,fm,q,cond)
  !This routine calculates the electrical conductivity from the perturbed distribution function.
  use precision,only:p_
  use constants,only: m,pi
  implicit none
  real(p_),intent(in):: dv,v(m),fm(m),q(m)
  real(p_),intent(out):: cond
  integer:: j
  real(p_):: sum

  sum=0._p_
  do j=2,m
     sum=sum+v(j)**3*fm(j)*q(j)*dv
  enddo
  cond=4*pi/3*sum
end subroutine get_conductivity


